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The Big Bang Observer (BBO) is a proposed space-based gravitational-wave (GW) mission designed pri- 
marily to search for an inflation-generated GW background in the frequency range ~ 10~^ Hz — IHz. The 
major astrophysical foreground in this range is gravitational radiation from inspiralling compact binaries. This 
foreground is expected to be much larger than the inflation-generated background, so to accomplish its main 
goal, BBO must be sensitive enough to identify and subtract out practically all such binaries in the observ- 
able universe. It is somewhat subtle to decide whether BBO's current baseline design is sufficiently sensitive 
for this task, since, at least initially, the dominant noise source impeding identification of any one binary is 
confusion noise from all the others (rather than instrumental noise). Here we present a self-consistent scheme 
for deciding whether BBO's baseline design is indeed adequate for subtracting out the binary foreground. We 
conclude that the current baseline should be sufficient. However, if BBO's sensitivity were degraded by a fac- 
tor 2 from the current baseline, then its ability to detect an underlying primordial background would depend 
critically on the value of pth, the threshold signal-to-noise ratio marking the boundary between detectable and 
undetectable sources. If BBO's sensitivity were degraded by a factor 4 from the current baseline, it could not 
detect a primordial background below ficw ~ 10^^^. 

It is impossible to perfectly subtract out each of the binary inspiral waveforms, so an important question is 
how to deal with the "residual" errors in the post-subtraction data stream. We sketch a strategy of "projecting 
out" these residual errors, at the cost of some effective bandwidth. We also provide estimates of the sizes of 
various post-Newtonian effects in the inspiral waveforms that must be accounted for in the BBO analysis. 

PACS numbers: 04.25.Nx,04.30.Db,04.80.Nn,95.75.Wx,95.85.Sz 



I. INTRODUCTION 

The Big Bang Observer (BBO) is a proposed space-based gravitational wave (GW) mission designed to search for stochastic 
gravitational-wave background generated in the very early universe QjQ. The design goal is to be able to detect primordial 
GWs with energy density i^qwif) 10~^^ in the frequency band 10~' Hz < / < 1 Hz. Standard, slow-roll inflation predicts 

^Gwif) < 10-16-10-15 11. 

To achieve this sensitivity to a primordial GW background, it will first be necessary to subtract from the BBO data stream 
the GW foreground generated by ~ 10^ — 10^ neutron star-neutron star (NS-NS), neutron star-black hole (NS-BH), and black 
hole-black hole (BH-BH) binary mergers, out to z ~ 5. This foreground "noise" has an amplitude substantially greater than 
BBO's instrumental noise, which in turn is probably substantially greater than the amplitude of the sought-for primordial GWs. 
To achieve BBO's goal, the GWs from the merger foreground must be subtracted to a level well below that of the primordial 
background. This means fliat flie amplitude of the residual, post-subtraction foreground must be < lO-^'^ of the pre-subtraction 
level. 

Will it be possible for BBO data analysts to subtract out the binary merger foreground to this accuracy? This question is non- 
trivial to answer precisely because confusion noise from unresolved mergers can in principle dominate the BBO noise spectrum. 
To decide which mergers are unresolvable, one needs to know the full BBO noise curve, including the level of confusion noise 
from the unresolvable mergers. But to determine the level of confusion noise, of course one needs to know which mergers are 
unresolvable. Clearly, one needs somehow to solve both these problems simultaneously. 

The focus of our investigations will be on NS-NS mergers, since these are the most problematic for BBO. The less numerous 
BH-BH and BH-NS merger events will have higher signal-to-noise ratios and therefore should be easier to subtract. If we find 
the NS-NS mergers can be almost fully subtracted from the BBO data stream, then the same should be true for the BH-BH and 
BH-NS mergers. 

How, in practice, will almost all the NS-NS mergers be subtracted out? We imagine that something like the following iterative 
scheme could be used: begin by resolving and subtracting out the brightest merging binaries (i.e., those with highest signal- 
to-noise-ratio), then resolve and subtract the next brightest ones, etc - regularly updating all the parameters of the subtracted 
binaries, as one goes along, to give the best global fit. Each subtraction decreases the foreground confusion noise and so 
increases the distance out to which NS binaries can be resolved. Will such a scheme suffice for BBO? The aim of this paper is to 
answer that question without actually having to carry out the whole procedure. We develop a method for determining the likely 
efficacy of foreground subtraction in a self-consistent manner Our method is (very roughly) as follows. Imagine that BBO is 
suiTounded by a huge sphere out to some redshift z, such that NS-NS mergers inside the sphere (i.e., at redshifts less than z) can 
all be individually resolved and subtracted (using realistic computational power), while none of the sources outside the sphere is 
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resolvable. This redshift z marking the boundary of the resolvable sources is not known initially, so we start with a reasonable 
guess. We then calculate the confusion noise due to all NS-NS mergers (NSm) at redshifts greater than z, <S'^^'"'^^(/), which 
we add to the instrumental noise S'^"''*(/) to obtain the total noise: 

5r(/) = ^r*(/)+5f (1) 

One can use this total noise level, S^^{f), to improve one's estimate of z, and iterate this procedure until z converges. 

Actually, of course, the detectability of any particular NS-NS binary depends not just on its distance (or redshift), but also on 
11 = L ■ N, where L is the normal to the binary's orbital plane and N points along our line-of-sight. (The binary's detectability 
also depends, of course, on the other three angles describing the binary's orientation and position on the sky, but to a much lesser 
extent.) Our calculation does properly account for the /i-dependence of the binary's detectability; i.e, we take 2 to be a function 
of /i, not a single number. 

We stress that there are actually two different sorts of confusion noise associated with merging binaries: the full signals from 
unresolved binaries (mentioned above), and the small errors that inevitably occur when waveforms from resolved mergers are 
subtracted out of the data. In Sec. IV we propose a method for dealing with these residual errors, by projecting out the subspace 
in which these errors can lie, at the cost of some bandwidth. We also estimate that this fractional decrease in BBO's bandwidth 
is small enough that for our purpose (deciding whether an iterative subtraction scheme is feasible) it can be neglected. 

We remark that our calculation is quite similar in spirit to a recent analysis of WD-binary subtraction in LISA data analysis, 
by Cornish et al. |4], which appeared when our own work was already at an advanced stage. In both cases, the idea is to use the 
requirement of self-consistency to arrive at a unique estimate of the efficacy of foreground subtraction, without actually coding 
up the whole analysis pipeline and testing it on simulated data. 

We also remark that a recent paper by Buonanno et al. 1 5 ] estimates that supernova explosions could provide another important 
BBO foreground, via the GW memory effect, but only if the anisotropy of neutrino emission is quite high, on average. For the 
rest of this paper we will neglect the possibility of a large foreground from supernovae. 

The organization of the rest of this paper is as follows. In Sec. II we give a brief overview of the BBO mission, its design 
sensitivity, and the foreground produced by merging NS binaries. In Sec. Ill we briefly explain why the most distant NS-NS 
binaries are effectively a noise source when it comes to resolving more nearby ones. In Sec. IV we summarize our proposed 
strategy of dealing with any residual subtraction errors by projecting them out. In Sec. V we provide estimates regarding the im- 
portance of eccentricity, NS spin, and high-order post-Newtonian (PN) effects in correctly subtracting out the resolved mergers. 
Besides being important for any future implementation of a BBO analysis pipeline, this catalog of effects is useful in estimating 
the threshold signal-to-noise ratio (SNR) pth required to detect NS-NS mergers. In Sec. VI we take a first cut at estimating pth, 
which we assume will be set by the then-available computational power. Our equations for self-consistently determining the 
efficacy of foreground subtraction are developed in Sec. VII. We solve these equations for a variety of assumptions regarding 
the NS merger rate, the detection threshold pth, and BBO's instrumental noise level, and display the solutions in Sec. VIII. We 
summarize our conclusions in Sec. IX. The derivation of one of the equations in Sec. VII is relegated to Appendix A. 

We use units in which G = c = 1. Therefore, everything can be measured in the fundamental unit of seconds. However, for 
the sake of familiarity, we also sometimes express quantities in terms of yr, Mpc, or Af©, which are related to our fundamental 
unit by 1 yr = 3.1556 x lO'^s, 1 Mpc = 1.029 x and IMq = 4.926 x lO^^s. 

For concreteness, we assume the universe corresponds to a flat Friedmann-Robertson- Walker model, with the universe's matter 
and vacuum energy densities being given by fini = 0.33 and i^A = 0.67, respectively. Our fiducial value for the Hubble constant 
isiJo = TOkms^^Mpc-^ 

II. OVERVIEW OF BBO AND THE NS-BINARY BACKGROUND 

A. BBO 

BBO is essentially a follow-on mission to LISA, the planned Laser Interferometer Space Antenna 0|, but optimized to 
detect GWs generated by parametric amplification during inflation. (For a review of inflation-generated GWs, see Allen f3l 
and references therein.) In the LISA band, 10^^ Hz - 10^^ Hz, an inflation-generated signal with 51gw ^ 10~^^ would be 
completely covered up by the foreground produced by galactic and extra-galactic white-dwarf binaries. By contrast, BBO will 
have its best sensitivity in the range ^ 0.1 Hz - 1 Hz. This band avoids the GW foreground produced by all the white dwarf 
binaries in the universe, which cuts off at / < 0.2 Hz (where the most massive of the WD binaries merge). In the BBO 
band, the dominant foreground GW sources are inspiralling NS-NS, NS-BH, and BH-BH binaries. BBO's baseline design, and 
corresponding instrumental noise curve, have been set in large part by the requirement that one must be able to individually 
identify practically all such inspiral signals and subtract them out of the data. An initial rough estimate suggested that the 
baseline "specs" in Table I are adequate for this purpose our primary task in this paper is to examine that issue much more 
carefully. 
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The current BBO design calls for four constellations of three satellites each, all following heliocentric orbits at a distance of 1 
AU from the Sun (see Fig.[0. Each 3-satellite constellation can be thought of as a "short-armed LISA". Two of the constellations 
overlap to form a "Jewish star"; the other two are ahead and behind by 27r/3 radians, respectively. Briefly, the idea behind this 
orbital geometry is that i7Gw(./) will be measured by cross-correlating the outputs of the two overlapping constellations in 
the Jewish star (much as LIGO attempts to measure flcwif) by cross-correlating the outputs of the Livingston and Hanford 
interferometers B)- The other two constellations give BBO its angular resolution: Ad ^ lO^^(SNR)^^ radians. It is not clear 
whether this angular resolution is strictly necessary for the purpose of measuring ^owif)^ but it will be immensely useful for 
BBO's secondary goal - to identify, map, and accurately determine the physical parameters of practically all merging compact 
binaries in the observable universe. 



FIG. 1: The Big-Bang Observer (BBO) consists of four LISA-like triangular constellations orbiting the Sun at 1 AU. The GW background is 
measured by cross-correlating the outputs of the two overlapping constellations. 



^From the output of each 3-satellite constellation (i.e., each "mini-LISA"), using time-delay interferometry (TDI) one can 
synthesize data streams that are free of laser phase noise and optical bench noise L2j JCLXU- A particularly convenient set of 
TDI variables to work with is {A, E, T}; all the GW information registered by each mini-LISA is encoded in these variables, 
plus the noises in these 3 channels are uncorrected with each other (i.e., they are statistically independent). Then, for instance, 
it is straightforward to find, for any source, the particular combination of {A, E, T} that yields the optimum detection statistic, 
and so to determine LISA'S optimum sensitivity to that source Holl . 

For our purposes, however, the following simplified treatment is adequate. As is clear from Fig. 4 of Prince et al. fl^, for NS- 
NS inspirals, each mini-LISA's sensitivity (using the optimum combination of the A, E and T channels) is practically equivalent 
to the sensitivity of two synthetic Michelson detectors, represented by the TDI variables X and Y . For our purposes, then, we 
can regard BBO, which is made up of 4 mini-LISAs, as formally equivalent to 8 synthetic Michelson interferometers. 

To construct the instrumental noise curve, S^^^{f), of each of these synthetic Michelson's, we used Larson's on-line "Sensi- 
tivity curve generator" 11211 . plugging in the parameters appropriate to BBO, which are listed here in TableU 



Symbol 



Value 



Laser power 
Mirror diameter 
Optical efficiency 
Arm length 

Wavelength of laser light 
Acceleration noise 



P 
D 
e 
L 
A 



3- 10" 



300 W 
3.5m 
0.3 

5 ■ lO'^ m 
0.5 fim 



TABLE I: BBO parameters. 



The parameters we adopt as reference values here are taken from the BBO proposal jll]; these parameters do not necessarily 
represent the latest thoughts on the mission's design (which is a moving target), but do provide a convenient baseline for com- 
parison. (Reference 1 1 ] also lists parameters for less and more ambitious versions of the BBO mission, referred to as "BBO-lite" 
and "BBO-grand", respectively, but in this paper we concentrate on the intermediate version, or "standard BBO".) In using the 
on-line generator, we have specified that the high-frequency part of S'j"*'' is 4 times larger than the contribution from photon shot 
noise alone. This is the same choice made in Fig. 1 of the BBO proposal [ IJ, and is consistent with assumptions typically made 
in drawing the LISA noise curve. As is conventional in the LISA literature, we take Sh{f) to be the single-sided, sky-averaged 
noise spectrum for each synthetic Michelson. This BBO instrumental noise curve is shown in Fig.|2l 



B. The NS-NS merger rate and the associated foreground noise level 



In this subsection we estimate the magnitude of the GW foreground from all NS-NS mergers. We denote the NS-NS merger 
rate (per unit proper time, per unit co-moving volume) at redshift z by h{z). The present-day density uq of merger remnants is 
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Frequency [Hz] 



FIG. 2: Shows the amplitude of the instrumental noise, \J /S'j"^'(/), compared to the amplitude of the (pre-subtraction) NS binary foreground 
(plotted for tiq ~ 10^^ Mpc^''yr~^) and the sought-for cosmic GW background (plotted for r2Gw(/) ~ 10"^^). Clearly, to reveal a cosmic 
GW background at this level, the NS foreground must be subtracted off, with fractional residual of < 10^^'^. 



related to n(z) by 11 



where 



H{z) = Ho v/a„(l + z)-^ + r!A . (3) 
As is conventional, we define r2Gw(/) to be the universe's fractional energy in GWs, per logarithmic frequency interval: 

"GW(/j = -77] 7^1 (4) 

Pc d(ln /) 

where pc = 3i?Q/(87r) is the universe's current energy density. Then the GW energy density (in the BBO band) due to (the 
inspiral phase of) all NS-NS mergers is given by fl3\ 

9 



= 1.7 X IQ-^^h-^ 



5/3 / r \ 2/3 



M Y'y f 



1.22MqJ VI Hz 

no \ /((l + z)-i/3) 



,103Mpc-VV 0-80 
The term ((1 + z)^^/'^) in Eq. ^ is the merger-rate-weighted average of (1 + z)^^^^, given by 



(5) 



OO 

((l + z)-i/3)^l / dz- 
no J (1 



(6) 



(l + z)4/3i7(2)- 



What is the universe's NS-NS merger rate history, ri{z)7 It is convenient to regard n{z) as the product of two factors: 

n{z) = no • r{z) , (7) 
where fiQ is the merger rate today and r{z) encapsulates the rate's time-evolution. 



5 



For r(z), we adopt the following piece-wise linear fit to the rate evolution derived in \U 

{l + 2z z <1 
1(5 -z) l<z<5 (8) 
z > 5 

For this r{z) and our fiducial cosmological model, one has 

no-no-(2.3-10iVr), (9) 

and ((1 + z)^^/"^) — 0.82. (We note that, as stressed in L13J . the value of ((1 + z)^^/'^) is actually quite insensitive to one's 
choice of the function r(z), generally being in the range ^ 0.7 — 0.9.) 

The current NS-NS merger rate, ng, is also usefully regarded as the product of two factors: the current merger rate in the 
Milky Way and a factor that extrapolates from the MiUcy Way rate to the average rate in the universe. The NS-NS merger 
rate in the Milky Way has been estimated by several authors; it is still highly uncertain, but most estimates are in the range 
10~^ — 10^'' yr~^ llSirifilllTll . To extrapolate to the rest of the universe, Kalogera et al. flol estimate that one should multiply 
the Milky Way rate by 1.1 — 1.6 x lO^'^ • h^Q Mpc^'^. That factor is obtained by extrapolating from the B-band luminosity 
density of the universe, and it is only a little larger than the extrapolation factor derived by Phinney in . Given the large 
overall uncertainty, in this paper we will consider 3 possible rates: ho — 10^^, 10^^, and 10^^ Mpc~'^yr~^. 

How many NS-NS merger events AiVm enter the BBO band during some observation time Aro? Summing the contributions 
from all redshifts, the rate N = AN„^/ Atq is easily shown to be 

/•CO 1 

N= 47r[aori(z)]2n(z)^dz, (10) 



where (for our fiducial cosmology) 



aoniz) = ^ f __=_4f___ (11) 

(12) 



Ho Jo ^{l-nA){l + z'f + 
dn 1 1 1 



Hol + z ^(l-rjA)(l + z)3 + flA 
This yields 



AiV.„ = 3.0 . 10^ ('^"l ( . ) (13) 



3yr/ VlO-^Mpc~V"^ 

The time required for a NS-NS inspiral signal to sweep through the BBO band will typically be comparable to BBO's lifetime. 
More specifically, the time remaining until merger, from the moment the GW frequency sweeps through /, is given (to lowest 
post-Newtonian order) by 

where A4 = ^^/^N'p/^ is the so-called "chirp mass" of the binary. (Here M is the binary's total mass and /i is its reduced mass.) 
Therefore, for two 1.4Af0 NSs, / « 0.205 Hz, 0.136 Hz, and 0.112 Hz at one year, three years, and five years before merger, 
respectively. 

Figure|3]plots the number of observable mergers during 3 years that occur closer than (any given) redshift z. We see that only 
^ 15% of mergers occur closer to us than z — 1. 

The (single-sided, sky-averaged) noise spectral density associated with any given GW background is liTJ: 

= 2"5^y3^C!w(/) (15) 



or 



[/^r (/)] = 8.8 X 10-- . h^o ( (^) . (16) 



6 




1 2 3 4 5 

Redshift z 



FIG. 3: The total number of NS-NS mergers closer than redshift z. The results here are normalized to a 3-yr observation period and lio = 
10"'^Mpc"Vr"^- 

The effective noise from all NS-NS inspirals (before subtraction) is plotted in Fig. 2, alongside the noise level from the sought- 
for inflationary background and BBO's instrumental noise curve. Clearly, the NS-binary foreground has amplitude ~ 10^ times 
higher than the (hypothetical) inflationary background's, in the BBO band, and so it must be possible to reduce (by subtraction) 
the foreground amplitude by more than ~ lO'^'^ to reveal an underlying primordial background. 

Given our r{z) and fiducial cosmological model, it is also straightforward to determine what fraction of S'^^'"(/) is due to 
sources farther out than some given redshift z. The result is plotted in Fig. |3 For example, 64% of the foreground spectral 
density is due to sources at z < 1, and 99% is due to sources merging at z < 3.6. Thus, very roughly speaking, one must 
subtract out all NS-NS mergers up to z « 3.6 to reduce the foreground noise amplitude by one order of magnitude. Of course, 
that conclusion is too simplistic, since the redshift out to which any particular NS binary can be observed depends on that 
binary's orientation as well as its redshift; see Section VI below for a proper accounting of this dependence. 




1 2 3 4 5 

Redshift z 



FIG. 4: Figure plots ^'^'^ j vs. i.e., it plots the fractional contribution of NS-NS binaries beyond redshift z to the total NS-NS 
foreground noise. 



III. UNDERSTANDING CONFUSION NOISE: WHY NS-NS CHIRPS INTERFERE WITH EACH OTHER 

So far, we have computed a spectrum for the NS-NS inspiral foreground, but we have not yet explained in what sense this 
foreground represents a noise source for BBO. We do so in this section, showing how GW signals from different mergers 
"interfere with" and so obscure each other. In this paper we simply sketch the main results; full details will be provided 
elsewhere L191 . 
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A. Brief Review of Optimal Matched Filtering 



Typical NS-NS merger signals will have amplitudes roughly two orders of magnitude smaller than the amplitude of BBO's 
instrumental noise. In practice, therefore, (some version of) matched filtering will be required to dig these buried signals out 
of the noise. Hence we will begin by briefly reviewing optimal matched filtering, partly to fix notation. For a more complete 
discussion, see |20] or i2lll . 

The output of N detectors can be represented by the vector SA{t), A = 1,2, ...,N. It is often convenient to work with the 
Fourier transform of the signal; the convention we use is 

/oo 
e'-'f'sA{t)dt, (17) 
-oo 

The output sa(0 is the sum of gravitational wave signal /ia(^) plus instrumental noise nA{t). In this section we will assume 
that the instrumental noise is both stationary and Gaussian. 'Stationarity' essentially means that the different Fourier components 
fiAif) of the noise are uncorrected; thus we have 

1 



nA(/) MfT = ^S{f - /')^h,AB(/), (18) 

where an overline ' ' denotes the 'expectation value' and Sh^Asif) is referred to as the spectral density of the noise. [When 
N=l (i.e., when there is just a single detector), we will dispense with detector indices and just write s{f) and 5h(/).] For our 
problem, we can restrict attention to the case where noises in different detectors are uncorrected; then we have 

1 



UAif) Mf')* = 2^(/ - /')^h,A(/)<5AB . (19) 

Given stationarity, 'Gaussianity' implies that each Fourier component has Gaussian probability distribution. Under the as- 
sumptions of stationarity and Gaussianity, we obtain a natural inner product on the vector space of signals. Given two signals 
gA{t) and kA{t), we define (g | k) by 



-< 



gl(/)feA(/)d/ 

^h,A(/) ■ ^ ^ 



It also follows from Eqs. (I19> and i20\ that for any functions gA{t) and kA{t), the expectation value of (g|n)(k|n), for an 
ensemble of realizations of the detector noise nA{t), is just (g|k). 

In terms of this inner product, the probabiUty for the noise to have some realization no is just 

p(n = no) cx e"<"°l"°>/^ (21) 

Thus, if the actual incident waveform is h, the probability of measuring a signal s in the detector output is proportional to 
/"^ ^ Correspondingly, given a measured signal s, the gravitational waveform h that "best fits" the data is the one 
that minimizes the quantity (s — h | s — h) . 

For a given incident gravitational wave, different reaUzations of the noise will give rise to somewhat different best-fit pa- 
rameters. However, for large SNR, the best-fit parameters will have a Gaussian distribution centered on the correct values. 
Specifically, let A" be the "true" values of the physical parameters, and let A" + AA" be the best fit parameters in the pres- 
ence of some realization of the noise. Then for large SNR, the parameter-estimation errors AA" have the Gaussian probability 
distribution 

p(AA") = A^e-^r^/^AA-AA^^ (22) 
Here Tap is the so-called Fisher information matrix defined by 



- \ 



(23) 



and M = \J det(r/27r) is the appropriate normalization factor For large SNR, the variance-covariance matrix is given by 



AA^AA/' (r-i)"'3 + ©(SNR)-! _ (24) 
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In the above notation, optimal filtering for some gravitational- waveform h{t) simply amounts to taking the inner product of 
h{t) with the data stream s{t). Assuming s = n + h, then 

(s I h> = (n I h) + (h I h) (25) 

The first term on the rhs of Eq. (125 > has rms value (h | h)^/^, so the signal-to-noise of the detection will be approximately 
given by 

SNR[h] = ^^^ = (h|h)i/2. (26) 
rms (h|n) 

B. Overlapping NS-NS chirps as a source of self -confusion 

Now imagine that the detector output s{t) consists of instrumental noise n{t) plus the sum of some large number of merger 
signals (labelled by "i"): 



(For simplicity, here we will consider the case of a single detector, and so eliminate the index A; the generalization to multiple 
detectors is trivial.) 

As explained above, optimally filtering the data for any particular merger waveform hj (t) is equivalent to taking the inner 
product (s I hj), which we can write as the sum of three pieces: 

(s I h,> = (n I h,> + ^(h, I h,) + (h, I h,) . (28) 

For the signal to be detectable, the third term should be significantly larger than the rms values of the first and second terms. 
We now explain why the second term can be sizeable; i.e., why different chirp signals can have substantial overlaps llssil . To 
simplify this discussion, let us use a slightly simpler version of the inner product; define 

/oo /•oc 
nnHfw^ / 9{t)Ht)dt, (29) 
-oo J —CO 

where the second equality in Eq. ( I29> is just Parseval's theorem. (Clearly, this is just our usual inner product, but without the 
"re-weighting by 1/ ^h" in the frequency domain. For white noise, where S^if) — constant, ( | ) and ( | ) are equivalent, except 
for an overall constant.) 

We now want to estimate the values of (n | h^) and (h^ | hj) for any two binary inspiral waveforms hi{t) and hj{t). In the 
nearly-Newtonian regime of interest to BBO, these are simple chirp waveforms: 

h,{t) = A,{t) cos <i>,{t), (30) 
hj{t) = Aj{t)cos't>j{t) , (31) 



where 



= cos / 27r/,(i')di', (32) 



= cos J 27rfj{t')dt' , (33) 

and where Ai{t), Aj{t), fi{t) and fj{t) are all slowly varying (meaning their fractional change during one cycle is << 1), 
and fi (t) and fj (t) are monotonically increasing. Then, since the integrand is highly oscillatory, it is clear that the integral 
J hi{t)hj{t) dt will show substantial waveform overlap only if there is some instant to when the two signals have the same 
frequency: 

Mto)=fj{to). (34) 

I.e., if one considers the "track" of each signal in the f-t plane, then tp is the instant of time when the two tracks cross. Using the 
stationary phase approximation, it is straightforward to show that I19ll : 

(h,|h,)« iA,(io)A,(io)|5/ri/2cos[A$o±7r/4], (35) 
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where A'I'o = [^i{to) — ^jito)]^ = [fi{to) ^ /j (io)]> ™d where the sign in front of the 7r/4 in Eq. ( I35t is positive when 
6f > and negative when Sf < 0. 
We want to use this resuh to estimate 

(h^lE^O' ^^^^ 

ie., to sum the contributions from all binaries whose f-t tracks overlap the j'^ track. Since the phase differences A$o at different 
intersections will clearly be uncorrected, the contributions accumulate in a random-walk fashion; i.e., the square of the sum is 
approximately the sum of the squares of the individual terms. Also, as we show in the next subsection, a typical NS-NS "track" 
will intersect a very large number of tracks from other merging binaries, so we are in the realm of large-number statistics. Finally, 
while the magnitude of each squared-contribution scales like | (5/ 1 ^ ^ , the number of terms in the sum scales like the average value 
of \df\, since the larger the "relative velocities" of the tracks, the more crossings. The dependence of the sum on the typical size 
of \5f\ therefore ends up cancelling out, and one can show the following 1 19]. Let H{t) — J^i hi{t) be the entire foreground 
generated by NS-NS chirps, and let iJ's spectral density be SB_{f), normalized so that 

/•oo 

W{i)= / ^H(/)d/. (37) 
Jo 

Then the expectation value of {hj \ )^ is given by 



(h, I ^ h, )' - i y h]{t)Sn{f,{t)) dt . (38) 

But the same result holds for the mean-square overlap of hj{t) with stationary, Gaussian noise n{t): 

(hTR^ = ^ / ■ (i)^h (/, it)) dt (39) 

with n'^{t) — Jp°° Shif) df. I.e., the mean-square overlap of a single chirp hj{t) with the chirp foreground H{t) (excluding hj 
itself) is the same as the mean-square overlap of hj (t) with stationary, Gaussian noise having the same spectral density as H. 
(It is straightforward to generalize this result to inner products with non-trivial frequency-weighting 1 19].) It is for this reason 
that in Eq. ([T) we simply add together the spectral densities of the instrumental noise and the "confusion noise" from unresolved 
chirps. 



C. The number of overlapping inspiral tracks in tlie f-t plane 

We saw in the previous subsection that two chirp signals have substantial overlap only if their tracks in the f-t plane intersect. 
Here we consider the track from a typical NS-NS inspiral and estimate how many other inspiral tracks it crosses. 

Let p{f) be the probability density of merger signals in frequency space; i.e., at any instant, p{f)Af is the average number 
of NS-NS GW signals received near the Earth that are in the frequency range [/ — A//2, / + A//2]. Since the BBO mission 
lifetime is vastly shorter than the age of the universe, we can assume p{f) is time-independent, implying 

p(/)^ = const = iV, (40) 

where, again, N = AA'm/ Atq is the total rate of mergers in the observable universe (from all z). The GW frequency derivative 
/ is given by 

df/dt = [Mil + z)f' . (41) 



so clearly p{f) oc f^^^^^. 

Now consider any one track in the f — t plane, and examine it in the neighborhood of some frequency /. It is easy to see that 
the rms rate Vc at which it intersects neighboring tracks is 

rc=0.5p(/)A/ (42) 
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where A/ is the rms variation in frequency derivatives for sources with GW frequency /. The 0.5 factor in Ea. J42t arises 
because, for any two neighboring tracks at any instant, there is a 50% chance that they are approaching each other and a 50% 
chance that they are separating. 

Using Eq. ( 14 1> . we see that the rms relative "velocity" of nearby tracks is 

^ = V3^. (43) 

/ McS 

where we define McS = M{1 + z), and where AA^off is the rms variation in this quantity. Now the fractional variation in 
A4 itself, AA4/A4, is probably at least of order 0.1. However, from Fig.|3lwe see that this is small compared to the variation 
A(l + z)/{l + z), which is - 0.4. 
Thus Tc is roughly given by 

re = 0.5(5/3)p/^^ ^ \n (44) 

independent of the particular frequency /. That is, the rate at which any particular track crosses all other tracks is about one-third 
the total merger rate from all observable sources, independent of where one is on the track. Thus, for any one track over the last 
3 years of inspiral, one expects of order 10^ crossings. This amply justifies our use of large-number statistics in the previous 
subsection. 



IV. CONFUSION NOISE FROM IMPERFECTLY SUBTRACTED WAVEFORMS 



NS binaries limit BBO's sensitivity to a primordial background in two ways. First, there will be some binaries that are too 
weak (because of their distance and/or orientation) to be individually identified and subtracted, and these "unidentified binaries" 
clearly represent a source of "confusion noise." Second, even identified NS binaries will not be removed perfectly from the data 
stream; inevitably (due to the finite signal-to-noise of the observations) there are subtraction errors, which represent a second 
source of confusion noise. This section addresses the confusion noise that results from subtraction errors. First we will prove 
a simple theorem regarding the magnitude of subtraction errors. Then we will sketch a simple strategy for largely eliminating 
their impact on other analyses by projecting them out, at the cost of some bandwidth. We estimate that lost bandwidth for BBO, 
and conclude that the loss is small enough that in the rest of this paper we can safely neglect it. 

We believe the analysis and strategy we outline here will also be useful in similar contexts, especially in dealing with problems 
of confusion noise in LISA data. Here we provide only a sketch of the main ideas; more details will be provided in a forthcoming 
publication LI 9.1 . 



A. Subtraction errors due to noise 



We have argued that, before searching for a primordial GW background, one will want to first subtract from the data the best 
fit to each identified inspiralling compact binary. However, because of detector noise, the best-fit values of the binary parameters 
will differ from their true values, and so the best-fit waveforms will be somewhat in error. What is the typical size of the error? 
That is easy to calculate: Let h{t) be some gravitational waveform immersed in noisy data, and assume the waveform depends 
on iVp physical parameters A" (a = 1, • • • , iVp). Because of the noise, the best-fit parameter values A" will differ from the true 
parameter values by I20I1 

JA" = A"-A"w (r-i)"''(n|a^h), (45) 
and, correspondingly, the best-fit waveform h{t) will differ from the true one by 

(5h = h - h 

(46) 

Using Eas.( l23> . (I24> . and ( I46> . we can immediately estimate the norm-squared of this residual error To lowest order in SX", 
we have 



' (47) 
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Thus the size of {Sh | Sh) is independent of the signal strength, but increases linearly with the number of parameters that need 
to be fit for. 

Eq. ( 147 > estimates the weighted integral of \6h{f)\'^; it says nothing about rms size of \6h{f)\'^ at any particular frequency /. 
Now, one can always calculate \Sh{f)\^ using (to lowest order) 

\6Hf)\' = dJiU)dih*{f) (r-i)"^ , (48) 

but for back-of-the-envelope calculations, it is reasonable to simply turn Eq. (147 > into a point estimate for the relative error: 



mf)\ 

\Hf)\ 



"(5h 









1/2 



SNR 



(49) 



For BBO measurements of NS-NS binaries, iVp 11 (cf. Section VI), and for a typical source (i.e, for a source at z « 1.5, with 
/i = 0.5, where /i is the cosine of the angle between the line-of-sight and the normal to the binary's orbital plane), SNR w 140, 
so5/i//i- 2.4 X 10-2. 

Given the extreme accuracy with which foreground sources must be subtracted, at first glance this level of error seems un- 
acceptable. However it would be a mistake to regard 5h as a completely random, additive noise source in the data. For one 
thing, after the best-fit signal h{t) has been removed from the data stream, the amplitude of noise plus residual is smaller (on 
average) than that of the noise alone. To see this, consider again the case of data s{t) = n{t) + h{t), and assume that the 
observation time is T, and that the data has been band-limited to [—/max, /max]- Then it is easy to show that the noise alone has 
squared-magnitude: 



(n|n) =2/maxT, (50) 

which is just the number of data points, for data sampled at the Nyquist rate 2/max- Next consider the magnitude-squared of 
the post-subtraction data set, s — h, where again h is the waveform that best fits the data. A straightforward calculation shows 
that lil 



(s-h|s-h> = (n - (5h I n - (5h> (51) 

= 2/maxr - iVp . (52) 

I.e., fitting out waveform h causes the norm-squared of the data to decrease below what is expected from noise alone. This is 
easy to understand: the fitting procedure takes out not only the signal h, but also that part of the noise that "looks like" the 
difference between h and some other waveform, h, having slightly different physical parameters. Stated geometrically, if one 
considers the manifold of physical gravitational waveforms, embedded in the vector space of possible measured signals, one sees 
that any piece of the noise that is tangent to the waveform manifold (at the location of the true signal) gets fitted out. Indeed, one 
sees from Eq. ( 145 > that it is just this piece of the noise, lying in the tangent space to the waveform manifold, that is "responsible" 
for the parameter estimation errors SX" in the first place. In the next subsection we outline a strategy projecting out this error 
before one searches for an inflation-generated background. 

Note that nothing in the above arguments required the signal to emanate from a single physical source. E.g., if H{t) is the 
entire foreground signal coming from sources, 

i/(<) =^/i,(t), (53) 

1=1 

and if each hi (t) is described by p parameters, then the full parameter space is described by Np — p x parameters, and 



(5H I 5H) =px Ns (54) 

to lowest order in 1/SNR. The total SNR^ of the foreground H is just Ng times the average SNR^ of the individual sources, 
and A^p is of course directly proportional to N^, so the fractional error in subtracting the whole foreground is just the fractional 
error in subtracting a typical source: 

5H/H^5h/h. (55) 

For BBO measurements of the NS-NS foreground, we thus estimate 5H/H ~ 2.4 x 10^^. 

As a digression, we remark that because our foreground consists of a large number of overlapping sources, it should not be 
surprising if there are some near-degeneracies that make it practically impossible to determine some of the physical parameters of 
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some of the sources. (These are cases where the affect on H{t) of adjusting the parameters of one source can be almost perfectly 
cancelled by adjusting the parameters of another source.) We bring this up to make the point that such near-degeneracies do 
not necessarily imply any degradation in one's ability to subtract out the foreground. Indeed, in the case of very high SNR (per 
source), it implies the opposite: near degeneracies would imply that the residual (5H is somewhat smaller than estimated above. 
The reason is simple: a near degeneracy means that the effective dimensionality of the signal space (near the actual signal) is 
smaller than the number of parameters being used to describe it. I.e., one could find a new parametrization using a fewer number 
of variables, N^. Then a repetition of the above arguments would yield ((5H | SH) = < Np. For the BBO case, where SNR 
per source is « 140, it probably will require detailed simulations to determine whether subtraction errors are larger or smaller 
than indicated by the high-SNR result, Eq. We leave this question to future work. 

B. Projecting out residual subtraction errors 

In this subsection we propose one strategy for effectively cleaning the BBO data of subtraction errors, after the NS-NS binaries 
have been subtracted out. Using this strategy, we argue that the impact of subtraction residuals (arising from instrumental noise) 
becomes sufficiently small that they can be ignored in the rest of this paper. We do not argue that our strategy is the best one 
possible, but rather offer it as an "existence argument" that some such strategy is possible. The use of any alternative strategy 
that leads to the same conclusion would not affect the main results of this paper 

The basic observation behind our strategy is that the residual 5H{t) is mostly confined to a surface within the vector space of 
all signals: the tangent space to the waveform manifold at the best-fit point. The corresponding errors in the subtracted waveform 
can be expanded in a Taylor series: 

5H{t) = daH{t)5\°' + ^da,dpH{t)dX"6X'^ + ■■■ (56) 

where a, /3 — 1, ...,p x Ns- The first-order piece on the rhs is the linear combination of iVp ~ p x Ns wavefunctions (the 
daHit)), with unknown coefficients (determined by the noise). We propose projecting these directions out of the data stream. 
This is simple in principle. Consider the operator 

P = /_(r-i)"/3|a„H)(9^H|. (57) 

where for simplicity we use here standard bra-ket notation of quantum mechanics. It is trivial to verify that = P and that P 
destroys any wavefunction of the form daH{t)6X" . We propose acting on the data streams with P before searching them for an 
inflation-generated background. 

What fraction of the data have we thrown away, by using P? For a fiducial 3-yr BBO lifetime, with, say, ^ 3 x 10^ subtracted 
sources, each determined by ^ 11 parameters, A^p ^ 3 x 10^. Assuming a 2-Hz sampling rate (sufficient for capturing most of 
the signal), with ^ 10* s of data and 8 independent channels, the dimension S of the full data space is 5 ^ 1.5 x 10^. Thus the 
fraction of the data that is discarded is only Np/S ^ 2 x 10^^, which is a negligible loss. 

So far we have discussed projecting out the first-order piece of the subtraction error; i.e., the piece linear in the parameter 
estimation errors SX". What is the magnitude of the second-order subtraction errors (i.e., the ones quadratic in SX")7 This is 
clearly given by 



{Sm I Sm) = -{dad^n \ d-ydeU) SX^SXf^X^SX' . (58) 

but evaluating the rhs of Eq. ( I58t is beyond the scope of this paper, and so we content ourselves with a cruder estimate. The 
second-order errors clearly scale like the square of the first-order errors, so a very crude estimate is 5^H/H ^ {SH/H)^ ^ 
6 X 10^^. Of course, this estimate is properly multiplied by some pre-factor (which can only be obtained by calculating of the 
rhs of Eq.|58]). Depending on this pre-factor and the actual level of the NS-NS foreground, these second-order subtraction errors 
could be comparable in size to the sought-for inflationary background. If this is the case, we would advocate projecting out the 
second-order errors as well. The second-order errors are linear combinations of second derivatives dadpH{t). It is important to 
notice that such second derivatives vanish identically unless a and (3 are parameters describing the same binary. Thus the vast 
majority of such second derivatives vanish. For each binary, there are (11 x 12)/2 = 66 non-vanishing second derivatives, so 
projecting out the second-order piece of the subtraction errors would cost only ^ 1% of BBO's bandwidth. A crude estimate 
of the size of third-order subtraction errors is 6^H/H ^ (SH/H)^ ~ 10^^. Clearly, unless the missing pre-factor here is quite 
large (of order 100 or more), it should not be necessary to project these third-order errors out of the data. 
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V. CATALOG OF RELEVANT PHYSICAL PARAMETERS AND RELEVANT EFFECTS 
A. Subtraction errors due to inaccurate waveform templates 

In the previous section, we outlined a method for handling subtraction errors arising from instrumental noise. Another poten- 
tial source of subtraction error is inaccurate theoretical template waveforms. Provisionally, we will regard a physical parameter, 
effect, or post-Newtonian term as "relevant for BBO" if neglecting it would lead to relative errors in our theoretical inspiral 
waveforms of size 5h/h > 10"^ (since errors of that magnitude could dominate over the inflationary background). Since each 
inspiral waveform contains ~ 10^ cycles, knowing the waveforms to 6h/h > 10~^ requires calculating the waveform phase to 
roughly one part in 10^^! 

The post-Newtonian (PN) expansion is clearly the right tool for constructing the waveforms, since the PN expansion parameter 
M/r is small in the BBO band: 

^.5.5xlO-fMl±ilV7_^V^ (59) 
r \ 2.8Mq J Vo.SHzy ' ^ ' 

where / is the GW frequency. If one uses PN waveforms, the only reasons for theoretical error would be 1) failure to calculate 
post-Newtonian corrections to sufficiently high order in the PN expansion, or 2) failure to account for all relevant physical 
parameters (e.g., the spins of the NSs). 

This section provides an initial "scoping out" of the questions of which physical parameters are relevant, and which post- 
Newtonian order is sufficient. 



B. Orbital Eccentricity 

1. Typical eccentricities of binaries in the BBO band 



Here we consider the implications of small (but non-zero) eccentricity for the subtraction problem. We begin by estimating 
typical eccentricities of NS binaries when they are emitting GWs in the BBO band. 

It is well known that radiation reaction tends to circularize the orbits of nearly Newtonian binaries. For small eccentricity e, 
decreases with the orbital period P according to oc pi^/^ |22]. For arbitrary e, the mutual scaling is given by L22] : 



p2/3 



pl2/19 



CX 



(1 



121 
304^ 



(60) 



The two known NS-NS binaries that dominate current merger rate estimates are PSR 1913H-16 and PSR J0737-3039. Extrapo- 
lating from today's values of e and P for these two binaries, using Eq. ( I60> . we estimate that their eccentricities when they pass 

-19/9 .2 - - n / f ^-19/9 

°0737 

- / \ \ -19/9 

examples, we will provisionally assume that typical eccentricities are '-^ [10^9 _ 10^^] I ^ j'^^ j . However we will also 
consider the implications of a subpopulation of NS binaries with considerably larger eccentricity. 



through the BBO band will be e^g^g « 4.6 x 10 ^ ( o 3Hz ) ^0737 ~ ^-O x 10 9 ^_X_^ , Based on these two 



2. Ejfect of non-zero eccentricity on waveform phase 



The effect of small, non-zero eccentricity is to slightly increase the inspiral rate; to lowest nontrivial PN order and to first 
order in e^, the increase (derivable from |22]) is given by: 



df/dt = d//dt|e=0 



1 



157 



(61) 



In the stationary phase approximation, we can write the Fourier transform of the emitted waveform (omitting tensor indices) 
as |23l 

(62) 



where ". . ." stands for higher-order PN corrections, and where the phase ^'(/) can be written as 

*(/) = *0(/) + *c(/). 



(63) 
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Here ^'o (/) represents the zero-eccentricity phase evolution and has the following PN expansion: 



^o(/) = const + 2nft, + - [%ttM{1 + z)f) '''^ 

20 /743 o/, 

H \ V - IGny^'^ + 

9 V336 4My ^ ^ 



(64) 



with y = {nAI{l + z) f)^^^, while ^c(/) re pres ents the phase correction due to non-zero e^, and is given (again, to lowest 
nontrivial PN order and to first order in e^) by ll23ll : 



*e(/) 



7065 



[TrM{l + z)] 



-5/3 2 i-19/9 f-34/9 



(65) 



187136 

Here cq is the binary's eccentricity at the moment that the GW frequency (more specifically, the frequency of the dominant, 
n = 2 harmonic) sweeps through some fiducial frequency /q. (Note that, by Eq. ( I60t . the combination Bq/q^^^ is a constant, to 
lowest nontrivial order.) 

Plugging in fiducial values, we can re-express Eq. (165 > as 



*c(/) 



0.21 



-0.3 Hz 



10 



-8 



{l + z)M 



1.22M, 







-5/3 



/ 



0.3 Hz 



-34/9 



(66) 



Note the very steep fall-off of 5'c(/) with increasing /. This J-34/9 f^H-off is much steeper than for the other PN correction 
terms in Eq. ( I64t . so it seems quite unlikely that errors in fitting for cq could be "absorbed" into compensating errors in the 
other parameters. While ^'c(/) is negligible for frequencies above a few Hz, it is typically of size ^ 27r at / = 0.1 Hz. Clearly, 
then, orbital eccentricity is a relevant parameter that must be accounted for, both in subtracting out individual sources and in 
projecting out residual errors. From Eq. ( I66> . we can also estimate roughly how accurately BBO can measure the eccentricity of 
each binary; it should be possible to determine Cq gjj^ to within A(eg 3 jj^) ~ [10~^/SNR] ~ 10^^". 



3. Contribution ofn 



'GW 



Non-zero orbital eccentricity implies that even the quadrupole piece of the gravitational radiation is no longer purely sinu- 
soidal, but exhibits harmonics at all multiples nv of the orbital frequency f (for integers n > 1). Let £"„ be the gravitational 
luminosity due to the rt"' harmonic. For small e. En oc e'^"^^' , so in the range of interest for e, only and Ei could potentially 
be significant. While both E^ and Ei are c>c e^, it is easy to show that the n = 3 contribution to ^^qI^ dominates over the n = 1 
contribution. Therefore we concentrate here on the n = 3 harmonic. 

The imoE3/E2 is Q 

^6„2 



E3/E2 « (3/2)''e^ , (67) 



from which one easily derives 
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'3\V3\ 



2 

6 .3 .13/9 

2J [2) 
1.6- 10~^^ 



(4)^'gw(/) (68) 

\ / / 2 \ \ / J \ -13/9 

no \ f (eo.3Hz) \ f ^ 



103MpcV V 10-8 J V0.3HZ 

where (cq gjj^) is the average value (for all NS-NS mergers) of at f — 0.3 Hz. For our fiducial estimate of (cq 3Hz)' 'his is 
significantly below the sought-for level of inflation-generated GWs, and so the extra harmonics generated by non-zero e can be 
neglected. 

However, our estimate that (cq 3 jj^,) ^ 10"^ was based on the few known examples of close NS-NS binaries; what if there 
is a subpopulation of NS-NS binaries that merge with substantially larger eccentricity (e.g., due to the Kozai mechanism |25||)? 
The ratio of the n = 3 to the n = 2 piece of the waveform, is clearly of order e. Thus the n = 3 piece must 

be subtracted (or projected out) if e > 10~^. Fortunately, as the previous subsection makes clear, if Eq 3 ^ 10~^, then the 
waveform itself will inform us of this fact, via the phase evolution of the n — 2 piece. 
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Unfortunately, to subtract h^^^, one needs to know both e and the perihelion angle cu (at some fiducial instant or frequency), 
since the latter clearly determines the relative phase of the n = 3 and n — 2 pieces. How accurately can wp.s Hz be extracted 
from the data? Since uj is encoded only in the n ^ 2 harmonics, we estimate that Awq.shz ^ min{7r, {eq ^ SNR)^^}. 
Hence, while the ft,"^"^ piece is relevant for Gq > 10~^, it will be impossible to subtract it when Gq 3 ^ 10^^ (since cjo.s Hz 
will be undetermined). Fortunately, even in this case, h^^^ can simply be projected out of the data (in the manner described in 
Sec. IV.B) since all possible realizations of h^^^{t) lie in a two-dimensional vector space. To see this, note that if all parameters 
except wo.sHz were known, then one could express h"^^{t) in the form A3(t)cos[3($3(t) + wo.shz)], where A3{t) and ^3(1) 
are both known functions, and this can be expanded as cos[3ti;o.3Hz] x ^3(0cos[3$3(t)] — sin[3wo.3Hz] x ^3(^)sin[3$3(t)]. I.e., 
h'^^^{t) is just some linear combination of two known waveforms, with (unknown) coefficients cos[3aJo.3Hz] and sin[3ti;o.3Hz]- 



4. Summary of effects of orbital eccentricity 

Extrapolating from the known NS-NS binaries, we have estimated that typical eccentricities for NS-NS binaries radiating 
in the BBO band will be e < 10^'*. At this level, they would have a significant impact on the phase evolution of the n = 2 
harmonic, but the n — 3 and n = 1 pieces of the waveform would be negligibly small. In this case, when projecting out residual 
errors, one need not worry about the perihelion angle oj. On the other hand, if some subpopulation of NS-NS binaries has 
^0 3Hz ^ 10^'^, then this will be completely clear from the data itself. For these binaries, both 3 and cjo.shz ^"e relevant 
parameters, to be used boh in subtraction and in projecting out residual errors. Finally, there are cases when cjo.3Hz is relevant 
but impossible to determine. Fortunately, even in this case, /i"^'^ can simply be projected out, at very modest additional cost in 
bandwidth. 



C. Spin Effects 

We turn now to the effects of the NS spins. Currently there are five known NS-NS binaries in our galaxy that will merge in 
a Hubble time (four binaries in the disk and one in globular cluster M15). In only one system-PSR J0737-are the spin periods 
of both NSs known. For PSR J0737, Pa = 22.7 ms and Pb — 2.77 s. In the other four systems, the radio-emitting neutron 
star is also a fast rotator, with P ranging from 28.5 ms to 59.3 ms. The fast rotators all have low spindown rates and so appear 
to be recycled pulsars. ^From evolutionary considerations, one expects exactly one of the companions to be rapidly rotating 
(consistent with what we find for PSR J0737). We estimate the effect of the bodies' spins on the gravitational waveform, for this 
presumed-typical case where one NS is rotating relatively rapidly (P ^ 30 ms), while the other is slowly rotating (-P > 1 s). 



1. Precession of Orbital Plane 

If the NSs are spinning, then the orbital angular momentum vector L does not have fixed direction, but instead precesses 
around the binary's total angular momentum vector J, due to an effective L x S coupling. When either 1) the two masses 
are nearly equal, or 2) the spin of one NS is much greater than the other, then the lowest-order precessional dynamics take an 
especially simple form- so-called "simple precession" tldt . In fact, we expect both these conditions to be satisfied in most NS- 
NS binaries, since (as mentioned above), we expect only one to be rapidly rotating, and since in those binaries where both NS 
masses are accurately known, the masses are indeed nearly equal. Therefore we shall use the simple-precession approximation 
to estimate the magnitude of precessional effects on the waveform. 

Following Apostolatos et al. |26], let Al be the precession amplitude; i.e., the angle between J and L. While Al depends on 
the magnitude and direction of the spins, the precession period depends on neither (to a very good approximation). The total 
number of precessions, from the moment the GW frequency sweeps through / until merger, is (for Mi « M2): 

It is useful to define dimensionless spin parameters Xi by Xi = \Si\/Mf. The Xi related to the spin periods Pi by 



where the are the NS moments of inertia. Label the faster-rotating NS "1". Assuming xi >> X2, the precession amplitude is 
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simply 

Al w 2.3 X 10"'*(1 - cos^eiLs)^^^ 

/ XI \ ( Mil + z) \ f f y (71) 
VO.Ol/ V 2.8Mq J \O.SUzJ 

where f?LS is the angle between L and 5*1. If we ignored spin-orbit precession when subtracting out the NS inspiral waveforms, 
we would make relative errors Sh/h ^ Al- This is < 10"'^ for Pi > 10 ms, and so these errors would typically be benign. In 
any cases where Pi is significantly less than 10 ms, this will generally be clear from the data (from its influence on the orbital 
phase evolution) and these very-high-spin systems would presumably be treated as a "special class", requiring more parameters 
to fit them than typically necessary. 



2. Effect of spin-orbit and spin-spin terms on waveform phase 

We next consider the effect of the spin-orbit and spin-spin interactions on the waveform phase. Since we have considered the 
effects of orbital eccentricity and orbital-plane precession in previous subsections, we simplify the analysis here by assuming 
that the orbit is circular and that the orbital angular momentum vector L and the two spin vectors, 5*1 and 5*2, are all aligned. 
Then in a post-Newtonian expansion of the waveform phase the lowest order terms involving the spin-orbit and spin-spin 

interaction are ll27ll 

^ r . (72) 

where the terms /3 and 7 are explicitly given by 

113 . 25M2\ ,^^,2r 



113 25 Ml 



(73) 

{M2/M)\L-S2)X2 



and 



Assuming Pi ^ 30 ms and P2 1 s, this implies xi ^ 0.01 and X2 4 x 10"^, and then /3 ~ 0.04, while \a\ ^ 2.5 x 10"^. 
So plugging in fiducial values (with AIi = M2 — 1.4Mq), the spin-related phase terms are 

-1/3 (") 

*.(/)--*xio-(^)(54) (1 + .)-"= 

In summary, the spin-orbit term /3 is clearly relevant, while spin-spin term a is negligible for typical cases. Thus, while it 
takes 6 parameters to describe (initial conditions for) the two spin vectors Si and S2, for typical cases the spins' influence on 
the waveform can be adequately subsumed into a single parameter, (3. 



D. High-Order post-Newtonian Effects, neglecting spin 

To-date, the post-Newtonian equations governing the inspiral of (quasi-)circular-orbit binaries have been derived through 
p3-5jg order beyond the lowest-order, quadrupole-formula level i28il . Is that good enough for accurately subtracting out the 
merger waveforms from the BBO data, or are even higher-order treatments called for? In this subsection, we do a rough estimate 
that suggests that the P'^^N equations are sufficiently accurate for this purpose (or are at least very close). Since we have 



17 



considered the effects of spin and orbital eccentricity in previous subsections, for this subsection we will specialize to the case 
of nonspinning NSs in (quasi-)circular orbits. 

We return again to the stationary-phase approximation of the waveform 



h{f)^{M{i + z)f'ry^[i + 

and to the PN expansion of the phase 



(76) 



= const + 27r/ic + - (87rAl(l + z)f) 



-5/3 



^ 20 /743 3/2 

IH \ y - IGiry^'^ 

9 V336 4M J ^ ^ 



(77) 



Terms up through P'^-^N have already been calculated. We want to estimate the size of the P^N term in the series, which 

corresponds to a term of the form |(87r7W(l + z)f)-^l^ x [(C + D{fi/M) + E{fi/MY H for some coefficients 

C, D, E, - ■ ■ . The coefficient C could be derived from the results in 129,1 ; we have not done that calculation, but it is clear from 
1I29I1 that C is of order lO'^. It seems reasonable to assume that the sum C + D{^/M) + E{^/M)'^ + ■ ■ ■ is also ^ 10^. The rest 
of the P'^N term, |(87rAl (1 + z) f)-^'^ y'*, has magnitude 



4.06 X 10 



_e ( M{l + z) 
\ 2.8Mq 



/ 
IHz 
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and so the full term is of order 10"'^ at / = 1 Hz. 

Thus the P^N contribution is just at the border of being relevant. We suspect the full P*N term will have been calculated 
long before BBO flies, but even today one could generate a "poor man's" P^N waveform by simply omitting the terms involving 
D{^/M), E{ii/MY, etc., but including the term cx C, which we repeat is easily derivable from published results. Because 
/i/M K, 1/4, the omitted terms could easily be an order of magnitude smaller than the C-term, and so would be truly negligible. 

Therefore we believe that already, today, one could produce PN waveforms that are sufficiently accurate for BBO, or that 
are at least quite close. However we add that if this view turned out to be too optimistic- if it did prove difficult to generate 
sufficiently accurate waveforms, corresponding to realistic solutions of Einstein's equation-then there is also an obvious fall- 
back strategy: use an enlarged space of "phenomenological waveforms," such as those developed by Buonanno et al. llsoll . to 
identify and subtract out the inspirals. The family of phenomenological waveforms would depend on a few more parameters 
than the physical waveforms, so projecting out subtraction errors would cost somewhat more bandwidth, but the estimates in 
Sec. IV.B show that this cost would still likely be minimal. Therefore as long as some member of the phenomenological family 
lies quite close to each true waveform, meaning 6h/h < 10^^, the phenomenological family would suffice for the purposes of 
inspiral-waveform subtraction. 



VI. THE DETECTION THRESHOLD pth 

The GW strength (at the Earth) of any NS-NS binary is characterized by its signal-to-noise-squared, p^. By p^, we mean 
the matched-filtering SNR^ for the entire 4-constellation BBO network (whose output is 12 independent GW data streams, 8 of 
which have good sensitivity to NS binaries). We want to estimate the threshold value p^j^ required for the signal to be detectable. 
There are basically two sorts of considerations here. If one possessed infinite computing power, then this threshold value would 
be set just by the requirement that one has sufficient confidence in the detection (i.e., that the false alarm rate be sufficiently 
low). However in practice we expect the search sensitivity to be (severely) computationally limited, which implies a somewhat 
higher detection threshold. 



A. Lower bound on pth set by the number of effectively independent inspiral templates 

Let A't be the number of independent templates required to cover the parameter space of NS-NS inspiral waveforms ('inde- 
pendent' in the sense that they have only modest overlap with each other). Then for a given threshold value pth, the number of 
false alarms generated by this entire set is ~ A't erfc(pth/'\/2) ~ A^t (27r)~^/^(pth)~^ e~''tii/2 in practice, one would probably 
want this false alarm rate to be no greater than ~ 0.01. How large is A't for our problem? This has not yet been calculated, but 
because pth depends only logarithmically on iVt, a very rough estimate will suffice for our purposes. 
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Consider the parameter space of 'typical' inspiral waveforms, normalized by (h | h) = 1. These are effectively described by 
10 parameters: 

A" = (A\...,AN) 

(79) 

= [^0, InA^cff, In/ioff, /?, el, $o, 0, (j), 6l, 0l] ■ 

Here, to is the instant of time when the (n — 2 piece of the) GW frequency sweeps through some fiducial value /o (e.g., 
/o = 0.3Hz); A^cff = A^(l + z) = m(1 + z); (3 is the spin parameter defined in Eq. (I73> (and approximated here as a 
constant); Cp is the square of the orbital eccentricity at to; $o describes the orbital phase (the angle between the orbital separation 
vector f and some fixed vector in the orbital plane) at tp; 4>) give the position of the source on the sky; and <})]_) give the 
orientation of the binary's total angular momentum vector L (which precesses slightly, but which we can typically approximate 
as constant). We have omitted from this list the perihelion angle loq and 5 of the 6 parameters characterizing the two NS spin 
vectors, since we estimated in Sec. V.C that they typically have a negligible impact on the waveform. The luminosity-distance 
to the source, Dl, has been omitted since it affects only the waveform's overall normalization. 

Now imagine covering our N-dimensional manifold of waveforms with a hypercubic grid, such that the overlap of any wave- 
form on the manifold with the nearest gridpoint is > (1 — x), where x is a number that characterizes the fineness of our grid. 
The number of gridpoints A^t is then 13111 

Nt « {N/Sx)^/^ JVfdXK.. dA^ (80) 

where F is the determinant of the Fisher matrix F^^ = (SqIi | df^h) (again, subject to the constraint (h | h) = 1). In our case 
N — 10, and we adopt x ~ 0.5 as our fiducial grid spacing, so (A^/8x)^/^ w 100. We can obtain a rough estimate of the 
integral J -\/rdA^ . . . dA^ from estimates of the sizes of the diagonal elements of F, as follows. For each parameter Aq, let 
nx^ — (5a„ \ h~^dhld\o\, where is the range of integration for the a*'' parameter and \ h^^dhld\a \ is supposed to represent 
some 'typical' or 'rms' value of this quantity. Then 



y VfdAMA2...dAN< 



riAi . . . tian . (81) 



The rhs represents a rough upper limit to the integral because it ignores possible cancellations in the determinant coming from 
the off-diagonal terms. Based on a post-Newtonian expansion of the waveform, of the form shown above in W Dl we derive the 
following order-of-magnitude estimates for the different factors: 



- 10^ 71inA^ ^ 10^ Jllnp ~ lo^ 
^ 10^, ~ 10^, n$o ^ 10^, no ~ 10^, 
nn,, - IQi (82) 

where = ngn^, n^j = ng^n^-^, and where we have used 5(3 ^ 0.5 and Sel ^ lO"*". Using the above estimates, we 
find TVt < 10'^^. Allowing for cancellations from off-diagonal terms, it seems reasonable to assume is in the range ^ 
10^° - lO^*", implying pth > 12.5 - 13.5. That is, if matched filtering reveals a NS-NS inspiral with total SNR > 13, then one 
can be confident it is not simply a randomly generated peak. 

Now, one could complain that we have undercounted iVt by restricting to the parameter space of 'typical' signals, whereas 
among the 10^ — 10^ NS binaries that BBO will observe, there are probably some atypical ones; e.g., binaries in which both 
NSs are rapidly rotating. And these must also be identified and subtracted, for BBO to do its main job. This complaint has some 
merit, but we do not dwell on it here, since in any case we expect that in practice pth will be set not by the false alarm rate, but 
by computational limitations. We turn to these next. 



B. Limitations due to Finite Computing Power 



^From the estimates in the previous section, one readily concludes that straightforward matched filtering for all templates in 
the template bank will not be possible. The simplest implementation would require of order ^ 10^ A^t floating point operations 
(since each year-long template has ~ 3 x 10® data points, if sampled at ^ 10 Hz). A well known, FFT-based trick to efficiently 
search over all tg 13211 reduces this cost by a factor ~ nt(,/[3 In(lO^)] 10®, but would still require computation speeds of 
^ 1028±3 flops (operations per second). Extrapolation of Moore's law to the year 2025 suggests that perhaps ~ 10^^ flops will 
be readily available, which is 1 1 orders of magnitude too small for the job. 
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Therefore one will need to devise a suboptimal (but computationally practical) search algorithm, and live with the attendant 
loss in sensitivity. It is beyond the scope of this paper to design such an algorithm and evaluate its efficiency. Fortunately, 
though, the problem of searching for NS-NS binary signals in BBO data is closely analogous to the problem of searching for 
unknown GW pulsars in LIGO data, and the problem of devising efficient search algorithms for the latter has been studied in 
some detail [33. 341 . We will estimate the threshold sensitivity of BBO NS-binary searches based on this analogy, so we digress 
to describe optimized LIGO searches for unknown GW pulsars. 

By unknown GW pulsars, we mean rapidly rotating NSs whose sky location, amplitude and polarization, and gravitational- 
wave frequency (at any instant) and frequency derivatives are all unknown, and so must be searched over. I.e., the unknown 
parameters are the sky location (6, (f>), four parameters describing the amplitude, polarization, and overall phase of the waves 
(these can be usefully thought of as two complex amplitudes-one for each GW polarization), and the gravitational wave fre- 
quency and frequency derivatives at any instant: /, /, /, /, etc. The typical magnitude of frequency derivatives is assumed 
to be d"//df" ^ //t", where t is some characteristic timescale (basically the NS's spindown-age), but these derivatives are 
otherwise considered independent. 

For GW pulsars, we briefly describe the most efficient schemes that have been considered to-date, which are semi-coherent 
and hierarchical (i.e., multi-stage) searches; we refer to Cutler et al. [34] for more details. A "semi-coherent" search is one 
where short data stretches (say, a few days long) are all coherently searched, using some technique akin to matched filtering, 
and then the resulting powers from the different stretches are summed. The method is only "semi-coherent" because powers 
are added instead of complex amplitudes; i.e., information regarding the overall phase of the signal in different stretches is 
discarded. This allows one to use a much coarser grid on parameter space than would be required in a fully coherent search of 
the same data. The basic idea of multi-stage searches is as follows. In the first stage one searches, semi-coherently, through 
some fraction of the data (say, a month's worth), and identifies promising "candidates" in parameter space. One then follows up 
these candidates in the second stage, using a higher resolution on parameter space (a finer grid) and more data. This generates a 
second, sublist of candidates, which one then investigates with even higher resolution and yet more data, and so on. The idea is 
to reject unpromising regions in parameter space as quickly as possible, so as not to waste valuable computer resources on them. 
After A^s semi-coherent stages like this, any remaining candidates are verified using a final, fully coherent follow-up search in a 
very tiny region of parameter space. A priori, the best value for Ns is unclear; it was shown in Cutler et al. Ii34ll that for realistic 
GW pulsar searches, the gains from increasing saturate at iVg = 3 semi-coherent stages. 

The GW signal from a NS binary is practically the same as the signal from a low-frequency GW pulsar (except the binary's 
orbital frequency changes on a much shorter timescale than the spin-period of slowly rotating NSs). In both cases, the signal 
is essentially monochromatic at any instant, with a frequency that is slowly time-varying. In both cases there is an unknown 
sky position, two unknown complex amplitudes (equivalent to D, Oi^, and $o in the NS-binary case). The optimal statistic 
for searching over the two complex (four real) amplitudes, in both the GW-pulsar and NS-binary cases, is the F-statistic, which 
follows a chi-squared distribution with 4 degrees of freedom 1 35, 36]. (The distribution is the same no matter how many detectors 
are combined in the analysis; the 4 d.o.f. correspond to the 2 complex-or4 real -unknown amplitude parameters. The fact that 
BBO is composed of 4 LlSA-like constellations outputting 12 independent data streams does not affect this counting.) The 
biggest difference between the two sources is that for actual GW pulsars, the signal's intrinsic amplitude can be approximated as 
constant over the observation time, while in the NS-binary case, the GW amplitude grows significantly during the observation 
time. However we do not consider this difference as very important when comparing detection thresholds, especially because 
the search sensitivity is really set by the early stages, where the coherent integration times will be significantly shorter than one 
year 

The sensitivity of the GW-pulsar search is limited by the size of the parameter space one wishes to search; e.g., for an all-sky 
search, the size of the parameter space is set by the maximum frequency /max and the shortest spin-down age T,nin that one 
wishes to search over. We now try to choose a search-space that makes the LIGO GW pulsar search comparable in difficulty 
to the BBO NS-binary search. The pulsar parameters (/, /, /) are closely analogous to the NS-binary parameters {Ai. /i, 
which control the inspiral rate. Assuming a search up to frequency /max — 1000 Hz, and an observation time of Tq = 1 yr, 
we estimate n^rijirfy ~ (/max21))'^(2o/Tmin)^ ~ 3 x 10'^^(1 yr/rmin)^- Using the estimates from Eq. ( I82> . we find that 
njUfU-f ^ n\nMn\nfinp forrmin ~ 300 yr. 

Continuing our comparison of the LIGO/pulsar and BBO/binary searches, we note that because three of BBO's four con- 
stellations have separations of order 1 AU (« 500 s), the number of distinct patches on the sky that must be searched over is 
~ (47r)(27r x 0.3 Hz x 500 s)^ ^ 10^. In comparison, for GW pulsar searches, the number of distinct sky patches is set by 
the Earth's rotation about its axis, and is ~ 3 x lO"*, or roughly 300 times fewer (This counting assumes that the larger, but 
more slowly varying, Doppler shift due to Earth's motion around the Sun can be absorbed into the unknown pulsar spin-down 
parameters, which should be true for integration times shorter than a few months. This is good enough for our purposes, since 
the sensitivity of the search is really set at early stages, where only a month or two of data is examined.) On the other hand, 
assuming sampling at ~ 10 Hz for BBO and sampling at ^ 3 kHz for the LIGO network, a year-long GW pulsar template 
contains ~ 300 times as many points as a year-long BBO NS-binary template, so each coherent integration requires about 300 
times more floating point operations in the LIGO/pulsar case than in the BBO/binary case. 

Therefore we conclude that a LIGO/pulsar search for unknown NSs, over a parameter range set by (/max = 1000 Hz, Tmin = 



20 



300 yr), is comparable in difficulty, computationally, to the BBO/binary search. The code used by Cutler et al. fsT] to calculate 
the efficiencies of multi-stage GW pulsar searches was re-run for this parameter range, assuming an available computational 
power of 10^^ flops (and computation time of one year). For this parameter range and computational power, LIGO/pulsar search 
with 3 semi-coherent stages (plus a final, coherent follow-up) should be able to detect GW pulsars with p as small as 20 1 3§] 
(with false-dismissal rate = 10% and false-alarm rate = 1%). Therefore we estimate that BBO will also be able to detect and 
remove NS binaries with p > pth ~ 20 (or roughly 50% higher than the minimum pth ~ 13 required for detection confidence). 

However: as in the last subsection, one could complain that we have counted only the cost of searching for 'typical' binaries, 
whereas in practice most of the computational budget may be spent on searching for the few atypical ones. Also, we have 
assumed (reasonably, we think, but without justification) that the computational cost of identifying all the individual sources 
is greater than (or at least comparable to) the cost of finding the combined best fit. Also, the comparison was made for a 
single false-dismissal rate (10%), whereas we imagine that, in actual practice for the BBO analysis, one would want to do the 
BBO analysis in stages, with an ever-decreasing FD rate. Also, actual BBO searches may be plagued by many more outliers 
than would be present for the purely Gaussian noise that our sensitivity estimates were based on, and this would increase the 
threshold. For all these reasons, and because our method of estimating pth ~ 20 "by analogy" was so crude in the first place, we 
will investigate the efficacy of NS-binary subtraction for a range of detection thresholds: pth = 20, 30, or 40. 



VII. EQUATIONS CHARACTERIZING A SELF-CONSISTENT SUBTRACTION SCHEME 



Fix the values of the merger rate fiQ (which sets the overall magnitude of 5^^'") and the detection threshold p^j^. We want 
to calculate what fraction of the spectral density of the NS-binary foreground cannot he subtracted. For simplicity, we will 
assume that all NSs have mass 1.4AfQ Then our method for self-consistently determining proceeds by the following steps. 

Step 1: Adopt some initial "guess" value Fq. Based on this guess, we obtain a corresponding guess for the total noise level: 

(^G, /) = sr'if) + Fa • sj^^'^if) (83) 

Step 2: Based on this total noise level, we determine the redshift z{p,), out to which a NS-binary with orientation p, = L ■ N 
can be detected. This boundary z{p) (separating detectable and undetectable sources) is determined by the equation (derived in 
the Appendix): 



where the function f{p) gives the dependence of the squared waveform amplitude on p 

5 
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/(.) ^ , (^ + ^^)^ + ^^^ = + ep^ + 1). (85) 



J dp [(l + ;,2)2+4^2] 


Note we have normalized f{p) so that J{p)dp = 1. 

Step 4: We compute the fraction F^ of the NS-binary foreground that is due to sources more distant than z{p). Based on 
Eqs. (10) and (12) in Phinney this fraction is easily seen to be 



^ 

in terms of the integral 

oo 

^(--^-/'^- (i+.y/W) - '''' 

z 

where H{z) and r{z) = n{z)/ho are given explicitly in Eqs. (|3} and (|8}, respectively. 

So far, we have given an algorithm for computing F{Fo), i.e., for iteratively improving our initial guess Fq. An initial guess 
Fq leads to a self-consistent solution if F{Fq) — Fq. Clearly, we can short-cut the iterative procedure simply by looking for 
fixed points of this function. I.e., our last step is 

Step 5: Plot F{Fq), and look for fixed points, i.e., values Fq such that F{Fq) - Fq ^ 0. 

Our results are displayed in the next section. 
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VIII. RESULTS 



As motivated in previous sections, we calculate the efficacy of foreground subtraction for 3 different values of the present-day 
merger rate density, fiQ = {10~*, 10~^, 10~^} yr~^Mpc~^, and 3 values of the detection threshold, pth — {20, 30, 40}. This 
yields 9 different results for the self-consistent F representing the fraction of the foreground noise amplitude due to undetectable 
(and hence unsubtractable) NS binaries. We calculate these results both for the "standard BBO" design sensitivity, S^"™^*'{f), 
shown in Fig. 2, and for a less sensitive version having S^^^*' = 4 x S^'"^^^{f), i.e., with 2x higher instrumental noise ampli- 
tude. As a shorthand, we will refer to the latter as "standard/2" sensitivity. Our main results are presented in Sec. Vlll.A. In 
Sec.Vlll.B we gain insight into our results by exploring which binaries (i.e., which z and /i) are undetectable, for different no 
and pth- Finally, in Vlll.C, we consider the case of a larger foreground, fiQ = 10^^ yr^^Mpc^^; although this merger rate is 
unrealistically high, this case provides a rather interesting illustration of our general method. 



A. Efficacy of Background Subtraction for BBO with Standard and Standard/2 Sensitivity 



In Sec. VII we showed that self-consistent F values are fixed points of the function F{Fg), where Fg denotes a "guessed" 
value for this fraction. For standard BBO sensitivity, we find that the solution F is practically independent of no, for realistic 
merger rates. Specifically, we find 

F40 = 0.0015 , 



where our notation is that F20 is the solution F for pth = 20, assuming the standard BBO instrumental noise level, and similarly 
for F30 and F40. Therefore standard BBO is sensitive enough that the NS-NS foreground can be entirely (or almost entirely) 
subtracted, independent of the merger rate or detection threshold (for realistic values of those quantities). 

Next we consider BBO with "standard/2" sensitivity. We denote by the self-consistent solution for pth = 20 and stan- 
dard/2 sensitivity, and similarly for F^f^ and F^q. For this case, the results do generally depend on no (unlike for standard BBO). 
Our nine results for F', corresponding to the nine combinations of {uq, pth). given in Table II. To illustrate how these results 
are derived, in Fig. [S] we show the function F{Fq) — Fq for each uq, and for fixed pth ~ 30. The entries in the second row of 
Table II are just the the Fq values where the three curves in Fig.|5]pass through zero. 



no 

Pth 


^inst ^ ^st.inst 


10"* 


10-^ 


IQ-^ 


20 


0.0015 


0.0015 


0.0015 


30 


0.071 


0.077 


0.11 


40 


0.15 


0.17 


0.55 



TABLE II: Results for "standard/2" sensitivity. Table lists F' = (5'NSm,>z^g,NSm^i/2 different combinations (no, Pth)- F' is the 
amplitude of confusion noise from unsubtractable NS binaries, divided by the total foreground amplitude. 



None of the F' values is Table II is zero; which ones are sufficiently small that unsubtracted binaries would not significantly 
interfere with BBO's main goal? To answer this, in Tablelllllwe give the ratio [5^ ^'°'>"(/)/5^W(j)]i/2^ evaluated at / = IHz, 
for each combination (no, Pth)- Again, S'i^'^"''>^(/) = (F'^S^^^^f), while in Table III [5^^(/)]i/2 is the noise spectrum 
for a primordial background with flQwif) = 10"^^. Thus, ratios smaller than one indicate that the unsubtracted piece of the 
foreground is smaller than a primordial background with this energy density. We see that if pth — 20 (i.e., if the detection 
pipeline can uncover almost all NS binaries with total SNR = 20), then even with standard/2 sensitivity, BBO would still be 
able to detect a primordial background having flcwif) > 10~^^. 

However Tablelllllalso shows that if pth = 30 or 40, and instrumental sensitivity is standard/2, then BBO would be unable to 
detect primordial background of JIqw (/) ^ 10" (since it would be "covered up" by the unsubtractable part of the foreground). 

We point out that entries for the case (no = 10"^, pth — 20) in Tables II and III should not be taken too literally, since in that 
case our solution F' corresponds to less than one unsubtracted binary. (A single merging NS binary at z = 5, even with /i = 0, 
contributes ^ 10" to our local flow if), in the BBO band.) What this means, of course, is that our solution F' lies outside the 
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FIG. 5: Shows the function F{Fg) - Fq for three merger rates: no = {10 10 ^, 10 "^jyr ^Mpc ^. All curves are for "standard/2" 
sensitivity and detection threshold pth = 30. 

range of validity of our equations, whose derivation implicitly assumed that at least one source was undetectable. Just as clearly, 
our main conclusions are unaffected. The proper interpretation of the {tiq — 10~^, ptii — 20) entries is that, for these values, 
BBO with standard/2 sensitivity would likely detect every single NS-NS merger occurring on its past hght cone. 



no 

Pth 


^inst ^ ^st.inst 


10"* 


10-^ 


10"® 


20 


0.030 


0.10 


0.30 


30 


1.4 


4.9 


22 


40 


3.0 


11 


110 



TABLE III: Table of ratios [S'^^'"'>^(/) /Sf ^(/)]^/^ evaluated at / = 1 Hz, for BBO with standard/2 sensitivity. Here 5*^ ^(/) is from a 
primordial background with S7Gw(,f) = 10~^^. Ratios smaller than one indicate that the unsubtractable part of the NSm foreground noise is 
smaller than this primordial background level. The results here are equivalent to those in Table II. 



We also repeated the above analysis for BBO with only standard/4 sensitivity, i.e, with S^^*'{f) = 16 • S";^*' 
This noise level is clearly inadequate, since even for pth = 20 and a low merger rate, no — 10"^ yr^^Mpc""^, we find 
j^NSm,>.-^^)/^GW(^)]i/2 « 3.0 at / = 1 Hz, for QGwif) = IQ-''- 

B. Further analyses of the subtraction scheme 

Here we expand on the results of the previous subsection, to improve understanding. In Fig.|6lwe plot the SNR of NS binaries 
having /i = (i.e, those seen edge-on: the least detectable case) as a function of z, under three different assumptions. The 
lowest curve (solid line) assumes standard BBO instrumental noise and assumes that the foreground confusion noise is the full 
S^^'^if) (i.e., the level before any subtraction), with ho = 10 ^yr ^Mpc In this case, assuming pth = 30, all binaries out 
to z « 1.5 could be detected, even without first subtracting out the brightest sources. (And of course, the binaries with more 
favorable orientations could be detected even farther out.) In an iterative subtraction scheme, one would begin by subtracting 
out all the high-SNR sources, which would lower the total noise and allow one to "look deeper" in succeeding iterations. For 
standard BBO, this iterative scheme reaches the point where there are zero, or almost zero, unsubtracted sources, and then the 
total noise is just the instrumental noise. 

The SNR for this "instrumental noise only" case is shown in the upper (dot-dashed) curve in Fig. |6l ^From this curve one 
sees immediately that F = is indeed a self-consistent solution: even the sources with = at z = 5 are detectable. What 
Fig. |6] cannof show is whether F = is the only self-consistent solution, but the rest of our analysis shows that this is true 
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FIG. 6: Shows the SNR ratio of NS-NS mergers with fi — in Eq. <84> . The total noise level is different for each curve. The solid curve is 
for standard BBO instrumental noise plus confusion noise from all NS binaries. The dotted curve is for "standard/2" instrumental noise plus 
foreground corresponding to no = 10~^ yr~^Mpc~'^, with F30 = 0.077 The highest curve is for standard-BBO instrumental noise and zero 
foreground noise. The horizontal line just highlights SNR = 30. 

(again, for standard BBO sensitivity and pth < 30). This has the practical implication that our envisioned iterative subtraction 
procedure should not get "stuck" at some higher F value: it can keep going until all binaries have been removed. The situation 
for standard/2 sensitivity is different, as illustrated by the middle (dotted) curve, which corresponds to the case pth = 30 and 
fiQ — 10~^ yr^^Mpc^^. For this case F' ~ 0.077, so the unsubtractable foreground noise is small compared to the instrumental 
noise, and the SNRs are roughly half the standard-BBO values. But then the ^ = binaries can only be detected to z w 2.2. 

The distribution of unsubtractable binaries, for BBO with standard/2 sensitivity, is explored further in Fig.0 which shows the 
maximal redshift to which NS binaries can be detected, as a function of /i. The three curves are for our three detection thresholds: 
Pth = 20, 30, 40; all assume no — 10^^ yr^^Mpc^'^. For pth = 20 (solid curve), only a tiny corner of the (z, /i)-space contains 
sources too weak to be detected, and the number of sources occupying that corner would be of order one (for uq — 10^^). For 
Pth = 30 or 40, the "undetectable regions" are clearly much larger, and contain several percent (or more) of all sources. 

5 



4 

N 3 

tn 
dJ 

cr 2 
1 



"0 0.2 0.4 0.6 0.8 1 

Orientation |i 

FIG. 7: Graph displays the maximum distance to which NS binaries can be detected, as function of their orientation angle n = L - N,foi three 
different detection thresholds pth- Here the instrumental sensitivity is "standard/2" and the merger rate density is no = lO^^yr^^Mpc"^ for 
all three curves. 



C. Confusion noise from a very strong NSm foreground 

The results for F in Eq. ( I88> had basically no dependence on the merger rate r'lo, and the F' results Table II showed only 
weak dependence on ng, except at the highest values of no and pth- The reason for this is simple: for BBO to succeed, the 
unsubtracted foreground noise must be smaller than the primordial background. Therefore, for BBO even to be "in the right 
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ballpark", the unsubtracted foreground must be well below the instrumental noise level. In this regime, the SNR of any source 
is set almost entirely by S^^^*', and so is insensitive to uq. Our results are consistent with the fact that, even with sensitivity 
degraded by a factor 2, BBO would still be "in the ballpark" (albeit insufficient for high pth)- 

However the dependence of F on fiQ becomes greater as one increases the merger rate, i.e., as unsubtractable binaries come 
to represent a significant fraction of the total noise. Because such cases display the full utility of our self-consistent method, we 
here show results for an unrealistically high merger rate: no = 10^^ yr^^Mpc^^. Fig.|8lshows the function F{Fq) — Fg for 
this no, for standard BBO instrumental noise, and for our 3 values of pth- Interestingly, each curve now has two zeroes; i.e., each 
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FIG. 8: Plots F{Fg) — Fq for an unrealistically high merger rate density, no = 10 ^ yr ^Mpc ^. In contrast to cases with lower no, each 
curve now has two zeroes. However only the higher zero (larger F value) can be reached by an iterative subtraction scheme. 

case has two self-consistent solutions. A moment's thought, however, convinces one that the larger of the two solutions is the 
only one that is accessible by an iterative subtraction scheme. Such a scheme essentially starts at the right-most end of the curve 
and proceeds along it, moving to the left as sources are subtracted, until it reaches the first zero of F{Fq) — Fq. At that point, 
any undetected source is too deeply buried in the noise of the other undetected ones (plus the instrumental noise) to be identified. 
Otherwise stated: while you can self-consistently "be at" the lower- i^Q solution, the class of schemes we are considering cannot 
"bring you there" (and we suspect that no scheme can). Thus, we see that in a universe with no = 10^^ yr^^Mpc^'^, more than 
half the foreground noise would be unresolvable by standard BBO. 



IX. CONCLUSIONS AND FUTURE WORK 



We have calculated the efficacy of an iterative procedure for removing merging NS binary signals from the BBO data stream, 
as required to detect any underlying, inflation-generated GW background. Our calculation basically required as inputs: a) the 
BBO instrumental noise curve, b) an estimate of the extragalactic NS-NS merger rate (as a function of z), and c) an estimate of 
the inspiral SNR required for detection, with realistic computing power. We find that the current design sensitivity is sufficient 
to allow data analysts to subtract out the merger waveforms, for the entire range of reasonable merger rates. If BBO were 
less sensitive by a factor 2 (meaning a factor 4 higher in S^^^*"), then BBO's success would depend on having a rather good 
detection algorithm, capable of finding almost all sources whose total SNR exceeds ^ 20. If BBO were less sensitive by a factor 
4, unsubtractable sources would simply "cover up" any underlying primordial background with JIgvv ^ 10~^^ (or somewhat 
lower if NS-NS merger rates are at the low end of the predicted range). 

Our goal was to estimate the efficacy of an iterative subtraction procedure, without actually trying to implement it. Of course, 
simulations of this procedure, to confirm our calculation or reveal holes in the argument, would also be very interesting. In 
particular, it would be important to confirm that our proposed projection technique on the cleaned data stream does sufficiently 
de-contaminate it of residuals from imperfect subtractions of resolved binaries, as we have assumed in this paper. A careful 
simulation of the BBO data analysis process would also lead to a firmer estimate of the threshold SNR pth required for merger 
detection in practice, as a function of available computing power. 
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APPENDIX A: DERIVATION OF EQ. ^ 

In this Appendix we derive Eq. ( I84> . We begin by averaging over all angles, including fj. = L ■ N; we return to the fi- 
dependence near the end. 

Consider first a single synthetic Michelson data stream from a single LISA-like detector. Let the waveform at the detector be 
hij{t) — h+{t)efj + hx where and e!^j are "+" and "x" polarization tensors, respectively. The average matched-filter 

SNR^ for some source (where the average is over source-direction and polarization angle) is given by 

Jo '^'h(7) 

where, as throughout this paper, S^if) is the "sky-averaged" noise spectral density. Parseval's Theorem states that 

1 

\h+{f)\'df^- hl{t)dt, (A2) 



2 „ 

and similarly for hy so for a chirping signal with a slowly changing frequency f(t), it is clear that 

\h+if)f + \hAf)\' = lihlit) + hi{t))dt/df , (A3) 

where the overbar denotes time-averaging. 

For now, consider some GW source at low redshift (z << 1). Then the rate at which the source loses energy due to GW 
emission is 

E{t) ^ AttD^ (7r/V4) < h-lit) + hl{t) > (A4) 
where D is its distance, and where the averaging is over all directions from the source. Therefore we have 

< \h+if)\^ + \hAf)\^ >^^^^5;S7^dt/d/. (A5) 

The product E{dt/df) equals |di?/d/|. For a circular-orbit binary, the energy is approximately 

E « -l^M/r«-i/i(M7r/)2/3 (^^^ 

= -ix5/^(7r/)2/3 (A7) 



from which we obtain 

1 



Using this result along with Eqs. ( lAU and ( IA5t . we arrive at 

The generalization of Eq. iA9\ to arbitrary redshift is accomplished by the standard replacement fyf\ M M.{1 + z) and 
D Z^L, where I?l is the luminosity distance. The /^-dependence of the waveform's strength-i.e., the /(/i) factor in Eq. \%A\ - 
follows almost immediately from, e.g., Eqs.(2a-2b) of I2al . Finally, to arrive at Eq. ( I84> . we multiply the rhs of Eq. \A9\ by a 
factor of 8, to account for the fact that at low-to-mid frequencies BBO is approximately equivalent to 8 independent Michelson 
detectors, each with the same noise density Sh{f)- 
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